
    #################################################################
    ### Table 3 - R code
    
    T1 <- 92
    T2 <- 114
    T <- length(T1:T2)  # number of congresses
    year.odd <- seq(from=1971, to=2016, by=2)
    

    ### Import the individual-level data

    Micro_Data1 <- read.csv("Table3_Data.csv", header=T)
    Micro_Data <- na.omit(Micro_Data1)
 
    ### Get the Congress-level data
    Macro_Data <- read.csv("Table_2_3_Congress_level_data.csv", header=T)
    names(Macro_Data)

    library(mvtnorm)
    library(R2WinBUGS)
    
    ### Setup data

    Y  <- Micro_Data$EP.extreme 
    Party <- ifelse(Micro_Data$party==100, 1, 2)
    GOP <- ifelse(Micro_Data$party==200, 1, 0)
    Dem <- 1- GOP
    Southern.State <- c(40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 51, 54) # VA, AL, AR (Arkansas), FL, GA, LA, MS, NC, SC, TX, plus KY and TN
    South <- ifelse(is.element(Micro_Data$ICPSR_state, Southern.State), 1, 0)
    Southern.Dem <- ifelse(Micro_Data$party==100 & is.element(Micro_Data$ICPSR_state, Southern.State), 1, 0)
    Northern.Dem <- ifelse(Micro_Data$party==100 & !is.element(Micro_Data$ICPSR_state, Southern.State), 1, 0)
    Southern.GOP <- ifelse(Micro_Data$party==200 & is.element(Micro_Data$ICPSR_state, Southern.State), 1, 0)
    Northern.GOP <- ifelse(Micro_Data$party==200 & !is.element(Micro_Data$ICPSR_state, Southern.State), 1, 0)
    X_energy.1 <- (Micro_Data$coal + Micro_Data$crude_oil + Micro_Data$natural_gas)/Micro_Data$population
    X_energy  <- as.vector((X_energy.1 - mean(X_energy.1, na.rm=T))/(sd(X_energy.1, na.rm=T)))


    X1 <- Dem
    X2 <- Micro_Data$Interior_cmt
    X3 <- Micro_Data$SEPC_cmt
    X4 <- South
    X5 <- X_energy
    X6.1 <- Micro_Data$Environment.extreme
    X6  <- as.vector((X6.1 - mean(X6.1, na.rm=T))/(sd(X6.1, na.rm=T)))

    # Congress-level IVs
    Divided <- ifelse(Macro_Data$President_party==Macro_Data$Senate_Majority, 0, 1)
    Majority_Margin <- Macro_Data$Senate_Majority_Seats/Macro_Data$Senate_Total_Seats
    gao_reports <- Macro_Data$gao_reports
    Oil <- log(Macro_Data$oil_real) 
    Warming <- Macro_Data$global_warming
    
    # Standardize the continuous variables
    Majority_Margin.std  <- as.vector((Majority_Margin - mean(Majority_Margin, na.rm=T))/(sd(Majority_Margin, na.rm=T)))
    gao_reports.std  <- as.vector((gao_reports - mean(gao_reports, na.rm=T))/(sd(gao_reports, na.rm=T)))  # GAO reports are already lagged.
    Oil.std  <- as.vector((Oil - mean(Oil, na.rm=T))/(sd(Oil, na.rm=T)))
    Warming.std  <- as.vector((Warming - mean(Warming, na.rm=T))/(sd(Warming, na.rm=T)))
    Warming.std.lag  <- c(Warming.std[1], Warming.std[-length(Warming.std)]) 
    
    Z1 <- Divided             
    Z2 <- Majority_Margin.std 
    Z3 <- Oil    
    Z4 <- Warming.std.lag
    Z5 <- gao_reports.std

    congress <- Micro_Data$time

    N <- length(Y)
    J <- length(Z1)

    n.beta <- 9
    n.gamma <- 6
   
    data <- list("N", "J", "n.beta","n.gamma", "Y", 
                 "X1", "X2", "X3", "X4", "X5","X6", "congress", 
                "Z1","Z2","Z3","Z4","Z5")

    # setup initial values
    init1 <- list(beta=rnorm(n.beta), 
                  beta2=rnorm(J), sigma.y=runif(1),
                  gamma=rnorm(n.gamma),
                  alpha=rnorm(J), sigma.alpha=runif(1))
    init2 <- list(beta=rnorm(n.beta),  
                  beta2=rnorm(J), sigma.y=runif(1),
                  gamma=rnorm(n.gamma),
                  alpha=rnorm(J), sigma.alpha=runif(1))
    init3 <- list(beta=rnorm(n.beta), 
                  beta2=rnorm(J), sigma.y=runif(1),
                  gamma=rnorm(n.gamma),
                  alpha=rnorm(J), sigma.alpha=runif(1))

    inits <- list(init1, init2, init3)
    
    # setup parameters
    parameters = c("beta", "beta2", "sigma.y", 
                    "gamma", "alpha", "sigma.alpha")
    
    Model.fit = bugs(data, inits, parameters, "Table3_bugs.bug",
                        working.directory=getwd(), bugs.directory="C:\\WinBUGS14",
                        n.chains=3, n.thin=5, n.burnin=2000, n.iter=12000)
    
    ### 4. Analyze the draws from the Posterior distribution ############################################

    library(coda)    
    chain1 <- read.coda(output.file="coda1.txt", index.file="codaIndex.txt", quiet=T)
    chain2 <- read.coda(output.file="coda2.txt", index.file="codaIndex.txt", quiet=T)
    chain3 <- read.coda(output.file="coda3.txt", index.file="codaIndex.txt", quiet=T)

    MCMC <- rbind(chain1, chain2, chain3) # Choose the chain that puts Democrats on the left

    iter <- nrow(MCMC)
    
    alpha.finder  <- seq(from=1, to=T, by=1)
    beta.finder  <- (T+1):(T+n.beta)
    beta2.finder <- (T+n.beta+1):(T+n.beta+T)
    gamma.finder  <- seq(from=(T+n.beta+T+2), to=(T+n.beta+T+1+n.gamma), by=1)
    
    n.beta2 <- length(beta2.finder)

    alpha.mc  <- MCMC[,alpha.finder]
    beta.mc   <- MCMC[,beta.finder]
    beta2.mc  <- MCMC[,beta2.finder]
    gamma.mc  <- MCMC[,gamma.finder]

    alpha.label <- paste(T1:T2, "th Congress", sep="")

    beta.label <- c("Democrats",
                    "Interior",
                    "SEPC",
                    "South",
                    "Fossil Fuel",
                    "Interior X Dem",
                    "SEPC X Dem",
                    "South X Dem",
                    "Fossil Fuel X Dem")
                     
    beta2.label <- year.odd

    gamma.label <- c("Intercept                                       ", 
                     "Divided Government",            
                     "Majority Seat Margin",
                     "Oil Price (log)", 
                     "Global Warming",
                     "GAO_reports"
                     )


    ### (1) Histograms of marginal posteriors of selected parameters
    ################################################################

    postscript("density_beta_trace.eps", paper='letter', horizontal=F)
    par(mfrow=c(3,2), mar=c(2, 5, 3, 3)) # 'mar -> c(bottom, left, top, right)'
    for(i in 1:length(beta.finder)){
    plot(density(beta.mc[,i]), main="", xlab="")
    title(main=paste("Marginal Posterior:", beta.label[i]), font.main=1)
    plot(beta.mc[,i], type="l", xlab="iteration", ylab="")
    title(main="Traceplot", font.main=1)
    }
    dev.off()

    ### Credible Intervals

    beta.CI  <- round(apply(beta.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    beta2.CI  <- round(apply(beta2.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    gamma.CI  <- round(apply(gamma.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    alpha.CI  <- round(apply(alpha.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    
    beta.all.CI <- beta.CI
    n.beta.all <- ncol(beta.all.CI)

    ### Plot the results

    label_beta <- beta.label
    
    rr <- 0.3 # adjust the x-axis of the label
    tt <- 0.5 # adjust the height of each panel
    x_lim1 <- -1 
    x_lim2 <-  1.5  #c(bottom, left, top, right)

    # use layout (instead of mfrow) to put three pictures in one figure
    postscript("Table3.eps", height=7, width=9, horizontal=F)
    par(mar = c(0,0,0,0), oma = c(3,4,0,0), mgp=c(1.5, 0, 0), xpd=NA, tck=-0.007)  #c(bottom, left, top, right)
    layout(rbind(c(1, 3), c(2, 3))) # this will put fig 1 on the top left, fig 2 on bottom left, and fig 3 on colunum 2

    plot(1:(n.beta.all), 1:(n.beta.all), type='n', col='black', xlim=c(x_lim1, x_lim2), ylim=c(1-tt, n.beta.all + tt),
            yaxt="n", xlab="", ylab="",  xaxt="n")
    #axis(side=2, las=1, at=1:(n.beta.all), labels=label_beta, cex.axis=.7)
    text(    beta.all.CI[4,] - rr, 1:n.beta.all, labels=label_beta,  cex=0.8) # use srt for text rotation 90=degree
    points(beta.all.CI[3,1:(n.beta.all)], 1:(n.beta.all), type='p', pch=19, col='blue') # Plot the Medians
    segments(beta.all.CI[1,1:(n.beta.all)], 1:(n.beta.all), 
             beta.all.CI[2,1:(n.beta.all)], 1:(n.beta.all), col='red', lty=1) # Plot CI
    segments(beta.all.CI[4,1:(n.beta.all)], 1:(n.beta.all), 
             beta.all.CI[5,1:(n.beta.all)], 1:(n.beta.all), col='black', lwd=3) # Plot CI    
    #title(main=paste("Senator-Level Coefficients"),  font.main=1)
    abline(v=0, lty=2)

    plot(1:n.gamma, 1:n.gamma, type='n', col='black', xlim=c(x_lim1, x_lim2),ylim=c(1-tt, n.gamma + tt), 
            yaxt="n", xlab="", ylab="")
    #axis(side=2, las=1, at=1:n.gamma, labels=gamma.label, cex.axis=.7)
    text(    gamma.CI[4,] - rr, 1:n.gamma, labels=gamma.label,  cex=0.8) # use srt for text rotation 90=degree
    points(gamma.CI[3,1:n.gamma], 1:n.gamma, type='p', pch=19, col='blue') # Plot the Medians
    segments(gamma.CI[1,1:n.gamma], 1:n.gamma, 
             gamma.CI[2,1:n.gamma], 1:n.gamma, col='red', lty=1) # Plot CI
    segments(gamma.CI[4,1:n.gamma], 1:n.gamma, 
             gamma.CI[5,1:n.gamma], 1:n.gamma, col='black', lwd=3) # Plot CI
    #title(main=paste("Congress-Level Coefficients"),  font.main=1)
    abline(v=0, lty=2)
    
    plot(1:n.beta2, 1:n.beta2, type='n', col='black', xlim=c(-0.5, 2), 
            xlab="", ylab="", yaxt="n")
    text(    beta2.CI[5,1:n.beta2] + rr, 1:n.beta2, labels=beta2.label,  cex=.8, srt=0) # use srt for text rotation 90=degree
    points(  beta2.CI[3,1:n.beta2], 1:n.beta2, type='p', pch=19, col='blue') # Plot the Medians
    segments(beta2.CI[1,1:n.beta2], 1:n.beta2,  
             beta2.CI[2,1:n.beta2], 1:n.beta2, col='blue') # Plot CI
    segments(beta2.CI[4,1:n.beta2], 1:n.beta2,  
             beta2.CI[5,1:n.beta2], 1:n.beta2, col='black', lwd=3) # Plot CI
    title(main="",  font.main=1)
    points(  rep(0, n.beta2+1), 0:n.beta2, type='l', lty=2) # use this instead of abline 
    
    dev.off()


